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The stability of a sheared magnetic field is analyzed in two-dimensional magnetohydrodynamics 
with resistive and viscous dissipation. Using a multiple-scale analysis, it is shown that at large 
enough Reynolds numbers the basic state describing a motionless fluid and a layered magnetic field, 
becomes unstable with respect to large scale perturbations. The exact expressions for eddy-viscosity 
and eddy-resistivity are derived in the nearby of the critical point where the instability sets in. In 
this marginally unstable case the nonlinear phase of perturbation growth obeys to a Cahn-Hilliard- 
QQ ■ like dynamics characterized by coalescence of magnetic islands leading to a final new equilibrium 

' state. High resolution numerical simulations confirm quantitatively the predictions of multiscale 

0^ ■ analysis. 
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I. INTRODUCTION 



Magnetic reconnection in two-dimensional magnetohydrodynamics (MHD) is one of the most intriguing problems in 
plasma physics, originally motivated by the observation of sudden and rapid release of energy occurring for instance in 
solar flares and tokamak disruptions. These phenomena are characterized by a change in the topology of the magnetic 
field lines and the formation of island- like structures. Resistive instabilities are often addressed as responsible of 
magnetic reconnection Their development is a consequence of a resistive boundary layer formed nearby a 

^vq ' magnetic field neutral line, where the ideal constraint of frozen magnetic flux is relaxed. Their growth rates typically 
. depend on fractional powers of the dissipative coefficients (through the dimensionless inverse Lundquist number 
and Reynolds number). When nonlinear effects are taken into account, the growth of resistive instabilities slows 
down, eventually reaching saturation ||. Numerical studies Q have shown that the nonlinear phase of evolution is 
characterized by the generation of large scale magnetic structures, by means of a coalescence process of small magnetic 
flux patches. 

In this Paper we present a new kind of instability which sets in when the influence of dissipative terms cannot 
be restricted to a limited domain and the boundary layer approximation does not hold anymore. This instability 
can be viewed as the MHD counterpart of the large scale hydrodynamical instability which is known to develop in 
highly anisotropic flows ||. Large scale instabilities are formally associated to the existence of negative values of 
the eddy- viscosity and/or eddy-resistivity. The exact values of the turbulent dissipative coefficients can be explicitly 
derived from the knowledge of the equations of motion and of the basic equilibrium state, by means of multiple scale 
analysis. In the case of a sheared magnetic field embedded in a motionless fluid, a negative eddy-resistivity instability 
can develop for high enough Lundquist/Reynolds numbers, causing the formation of a chain of magnetic islands. As 
observed for resistive instabilities, also in this case, during the nonlinear evolution, the number of magnetic islands 
decreases, due to the merging of pairs of close islands. This coalescence process eventually leads to the formation of 
a single island, whose width is determined by the largest scale allowed in the system. 

In Sectio n O , we briefly introduce the fluid equations of magnetohydrodynamics and their basic equilibrium states. 
In Section fn\ , the behavior of large scale perturbations is investigated making use of multiple-scale analysis and 
exact expressions for eddy-viscosity and eddy-resistivity are derived. The main result is the appearance of large scale 
transverse instabilities, associated to negative values of the renormalizcd dissipative coefficients. In Section IV, we 
focus on the case of marginal instability for which we obtain the effective equation for the large scale behavior and we 
show that the full nonlinear regime is characterized by the evolution towards a fixed point. In Section^, we present 
the results of direct numerical simulations (DNS), which display a clear quantitative agreement with the predictions 
of multiple-scale analysis, both in the linear and in the nonlinear phase. 



II. MHD EQUATIONS AND BASIC EQUILIBRIA 

MHD equations are relevant in many different physical contests, such as astrophysics, laboratory plasma physics or 
magnetized liquid metal dynamics. The applicability of this model, which is a fluid description of plasma, relies on 
the assumption that all the lengthscales under consideration must largely exceed the ion Larmor radius. In a strong 
external magnetic field which is oriented along z, B z 3> B±, the motion becomes almost two-dimensional and the 
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MHD model is well approximated by the 2D MHD equations for the magnetic flux function tp associated to the planar 
magnetic field (Bj_ = e z x V-0) and for the stream function tp of the incompressible planar flow (vj_ = e z x \/(p) ||: 

^ + =?7(V 2 V- Jo), (1) 

^+^,vV] =^,vV]+^vV, (2) 

where, following a standard notation, the convective terms are written as Jacobian operators ([/, g] = d x fd y g — 
d x gd y f). The forcing Jo = V 2- 0o represents an input of magnetic energy preventing the decay due to resistive (rj) 
and viscous {v) dissipation. 

The above equations have been normalized with respect to the characteristic macroscopic length L, magnetic field 
B and Alfven time ta — L/va, while 77 and v are, respectively, the inverse Lundquist number (rjc 2 /4ttvaL) and the 
inverse Reynolds number (v/vaL). The Alfven velocity va = B /(47rmn) 1 / 2 , is the velocity of small amplitude waves, 
propagating along the magnetic field B in a uniform plasma with density nm. 

The equations (HB) admit a basic equilibrium state ip = "00 = J~(x), ip = 0, where T is any function. We consider 
T[x) = cos x, in order to investigate the stability of a sheared magnetic field in a motionless conducting fluid, in 
a slab geometry with periodic boundary conditions. Such a configuration is widely used to study the evolution of 
reconnecting modes. It is known j| that the slab geometry can be trusted as realistic and resistive instabilities can 
develop only if the slab aspect ratio (L x /L y ) is less than one. For this reason it is of great interest the question related 
to the stability of the above configuration with respect to transverse perturbations on scales larger than the typical 
scale of the basic magnetic flux. These large scale instabilities are eventually responsible for the onset of the inverse 
cascade of the square magnetic potential A = (ip 2 ) in fully developed 2D MHD turbulence f|. 



III. MULTIPLE-SCALE ANALYSIS 



As a first step, using a standard multiple-scale technique, we derive analytically the large-scale magnetic flux ^ 
and stream function $ equations. The multiscale method is based on the idea of exploiting the separation of scales 
as a perturbative parameter. We consider a periodic box width Lb ox much larger than the basic magnetic flux typical 
length-scale (Lbox 3> L x ) in order to follow the dynamics on spatial scales of order 0(l/e)L x . Beside the fast variables 
(x, y, t) on which the basic flow evolves, a set of slow variables (X — ex, Y — ey, T — e 2 t) can be introduced. According 
to this choice the differential operators appearing in (jl|,g) are transformed to 

% -> di + eV, , d t -> d t + e 2 d T . (3) 

Expanding perturbatively in e the fields we obtain: 

1& = V (0) (3, X, y,Y,t,T) + (x, X, y, Y, t, T) + e 2 ^ 2) ... 

ip = p<°) (x, X,y,Y, t, T) + stpW (3, X, y, Y, t, T) + e V 2) .. . (4) 

The scaling of the slow time T is suggested by physical hints: we are looking for a diffusive behavior of large scales 
which takes place on times 0(e~ 2 ). It is worth noticing that in general the large-scale MHD dynamics is first order 
in time and space (the well-known a effect) , but it can be shown that this is not the case for parity- invariant basic 
configurations ||. 

By substituting (|^) and (Q) in (Q) and (j|) and by equating the same power of e, one easily finds a hierarchy of 
equations in which perturbations belonging to different order of expansion appear coupled and depend on fast and 
slow variables. The dependence on the fast time variables can be discarded by observing that it reduces to a transient 
not affecting the long-time behavior, thanks to the fact that the forcing and the basic flux are time-independent (a 
rigorous proof needs the construction of a Poincare inequality) . We look for solutions with the same periodicities of 
the basic magnetic flux in L\, ox . At each step we can distinguish the pure large field contribution from the small scale 
oscillating part: 

tpW = ^f( fe ) (x, Y, T) + ij (k) (X, Y, T, x) , p {k) = $ (fe) (X, Y, T) + ip ik) (X, Y, T, x) . 

Equations have to be solved recursively because solutions of lower order appear as coefficients in the following steps 
of the hierarchy. At each order one has to test the validity of the solvability condition. Indeed, the equation for the 
large scale magnetic flux is obtained as solvability condition at order e 2 , while the equation for the large scale 
vorticity V 2 $(°) comes out at order e 4 : 
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Let us focus our attention on diffusive terms in (^|) and (||): as a consequence of the anisotropy of the basic 
small-scale flow the eddy-diffusivities are anisotropic too. For longitudinal perturbations ($(°) = T),^ ^ = 

ty(°\X, T)) both viscosity and resistivity are left unchanged. On the other hand, for transverse perturbations ($(°) = 
^,(0) T),\I/( ) = (Y, T)) the renormalization of resistivity due to the small scale magnetic energy holds a negative 
term (—l/2v), while molecular viscosity is increased by the eddy contribution (l/2rj). 

It is interesting to notice that the possibility of a negative eddy-resistivity due to small scale magnetic energy has 
already been presented in ||. That result, obtained in the framework of closure approaches to 2D MHD turbulence, 
was suggested as an explanation of the square magnetic potential inverse cascade. Analogous results are reported 
in |u|, where the effective resistivity is shown to become negative in a small scale turbulent plasma, as long as the 
magnetic energy exceeds the kinetic one. 

By inspection of equations (|^) and (||), we stress that the large scale magnetic flux and stream function are 
linearly decoupled. If we assume large scale perturbations of the type ~ exp(r^T + iKxX + iKyY) and 

<j>(°) r*j exp(r$T + iKxX + iKyY), stability analysis leads to the following dispersion relationships 
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The stability problem can be tackled more easily by introducing the parameters P = l/2i]v and T — (Kx/Ky) 2 - 
Marginal stability lines (T* = T$ = 0) are then given by 

1 -P + T = , 

I + P + 2(1 - ^P) T + T 2 = . 

and plotted in Figure [l| 

For high enough values of molecular resistivity and viscosity (P < 1), the basic flow is stable against any large scale 
perturbation. Increasing the Reynolds numbers, the first instability sets in at P — 1 and T = 0, i.e. for transverse 
perturbations. We notice that for 1 < P < 16/9 the large scale vorticity is always stable (r$ < 0) and the magnetic 
potential growth rate is maximum in the case T = 0. 

In the following, we investigate the linear and nonlinear evolution of large scale perturbations in the neighborhood 
of the critical point (P = 1,T = 0), where we expect that there is no amplification of kinetic energy but that due to 
non-linear effects of coupling between velocity and magnetic field. 



IV. MARGINAL INSTABILITY 



Special attention deserves the development of this large scale instability for Reynolds numbers close to the marginal 
stability threshold. In this regime, it is possible indeed to follow the full nonlinear evolution of the perturbation. We 
show that the instability eventually reaches a fixed point characterized by a magnetic island of the size of the box. The 
nonlinear evolution in the marginal regime is described by two coupled equations which generalize the Cahn-Hilliard 
equation, found for the hydrodynamical counterpart of this system, the so-called Kolmogorov flow [p| Jll| , [l2]| . 

Let us suppose to move the parameters just above the marginal stability line: 

V = Vc (l - e 2 ) , ^ = ^ c (l-e 2 ) (9) 

where r\ c v c = 1/2 (P = 1). The perturbative parameter e is thus fixed by the distance between 77, v and their critical 
values rj c ,v c . We will take into account only transverse perturbations, since, as shown in Figure [TJ, the large scale 
magnetic flux linear instability is mainly transverse, close to P — 1. According to M), the transverse eddy-resistivity 



3 



in the neighborhood of the critical line defined by ([)]) is of order 0(e 2 ), thus suggesting a scaling for the slow time 
T — e i t. The decomposition rules (0) become 



d x -> d x d y ^ d y + ed Y , d t -» d t + e d T ■ 

The same multiscale technique described above can be adopted to solve perturbatively (|l|) and (||). 
At first order in e we obtain: 

ip = cos x + * (0) (Y, T) + e* (1) (Y, T) 

<y5 = 2?7 e £ sin a: . 

or 



(10) 

(11) 
(12) 



We notice that, according to the conclusions drawn above, the large scale stream function is linearly stable and it is 
simply driven by the magnetic flux. 

The evolution equations for and emerge as solvability conditions at order e 4 and e 5 : 
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5 T $(o) = - Vc d 4Y ^ i0) - 2r lc d Y Y^ {0) + 1277 c a YF * (0) (9y* (0) ) 2 
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(13) 
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In the first equation ( |l3| ) , one easily recognizes the renowned Cahn-Hilliard equation |llj] , which may be written in 
variational form: 



dt 



The existence of the Lyapunov functional 
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(15) 



indicates that asymptotically the solution of ([T^) in a bounded domain reaches a fixed point. This stationary solution is 
approached by a sequence of metastable states of decreasing dominating mode. We thus expect to observe a nonlinear 
evolution dominated by a magnetic island coalescence, analogous to the vortex pairing in 2D hydrodynamics |T^| . 

Equation ( p~4| ) is linear in 'pW, with coefficients depending nonlinearly on It also can be written as a gradient 

flow, with a Lyapunov functional 
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We conclude this section by observing that the dispersion relation for *(°) - exp(IT + iKY) now reads: 

27 

r = -—r, c K 4 + 2r, c K 2 . 



(16) 



(17) 



It implies instability (r < 0) for any (K < 0.77). We notice that information about the characteristic scale of unstable 
modes were totally absent in the general treatment presented in the previous section (see equation (|7]). 



V. NUMERICAL RESULTS 



The results obtained in the previous section have been checked by extensive direct numerical simulations of MHD 
equations In order to force a transverse perturbation, we integrate the equations on a rectangular slab with 

Lbox — L x — 2tt and L y > L x . In this way large scale instability can only develop on the y direction for an aspect 
ratio r — L x /L y < 1. 

Given the numerical values of parameters r\ and v, from (H) and the condition r\ c v c = 1/2 we have 
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b < 18 > 

which are used for the theoretical predictions of the previous section. 

The simplest check of our predictions concerns the growth rates of the instability which, in the initial linear regime, 
are given by the dispersion relation (|l7|). In physical (not rescaled) variables ( |l7| ) becomes 

27 

7 = -— Vc k 4 + 2 Vc e 2 k 2 (19) 

which shows that the largest unstable wavenumber is k max ~ 0.77e. The smallest transverse wavenumber is ki = r 
thus, in order to numerically observe the instability, it must be r < 0.77e. 

In Figure § we report the growth rates of the first modes for a simulation with r = 1/64, v = 0.49 and r\ = 1. 
We have e ~ 0.1 and thus only the first 4 modes are unstable. The initial perturbation is small, random and on 
all the first 20 modes, thus we are able to observe also negative 7's (stable modes). The comparison with the linear 
prediction ( |l9| ) is very good even for e not very small. The numerical data of Figure ^| are obtained by a linear fit of 
the logarithm of the mode amplitude versus time in the early stages of the simulation. 

Let us now consider the nonlinear stage of the perturbation growth. We describe here a different simulation with 
r = 1/16 and e ~ 0.32 which was advanced for a very long lapse of time. The nonlinear evolution will ultimately lead 
to a fixed point by a succession of long lasting quasi equilibrium states of decreasing wavenumber. The evolution of 
the amplitudes of the first 5 transverse modes computed from the direct numerical simulation is plotted in Figure ||. 
Observe that in this case the fifth mode k = 5/16 is linearly stable but it is non-linearly driven by smaller wavenumber. 

The typical linear time is now rather short, I/7 ~ 0(1), and the final stationary state, dominated by the largest 
available mode k\, is reached at very long times, t > 1000. At intermediate times, almost stationary metastable 
states, characterized by decreasing leading mode, are punctuated by fast coalescence processes. Most of the energy 
dissipation takes place during this fast reconnection processes. In Figure ||we display the period-two metastable state 
at time t = 200 and the final, period-one state at t = 20000. The dynamical picture arising from Figures || and ^ 
qualitatively agrees with the dynamics described by the Cahn-Hillard equation |1^] . 

To check quantitatively the validity of non linear multiscale analysis we have numerically integrated the Cahn- 
Hilliard equation for the large scale magnetic flux (|l^) with the same parameters of the DNS. As shown in Figure 
^| we find an impressive agreement even for very long times. The final relative amplitudes of the most energetic 
transverse modes is recovered within a 10 % accuracy. 

As a further test of the multiscale predictions, we checked the relations (JXl|-p~2|) during the evolution. At leading 
order in e, ^^°\y,t) is obtained by subtracting the basic flow cosx from the magnetic flux ip(x,y,t). The resulting 
field, which reveals to be indeed x- independent, is then used to reconstruct the stream function by means of (p^). 
The results for the configuration of Figure [| is shown in Figure ||. 



VI. CONCLUSIONS 



We have investigated the issue of stability of highly anisotropic, magnetic-energy dominated equilibrium states of 
the MHD fluid model equations. These configurations are known to be unstable for small values of resistivity leading 
to the formation of thin boundary layers in the nearby of the neutral line of the magnetic field. At variance with the 
above case, we focused our attention on the range of moderate Lundquist/Reynolds numbers, where the boundary 
layer approximation is not fruitfully applicable. In this situation we have shown analytically, by means of multiple- 
scale analysis, that large scale instabilities can arise above a definite threshold, and that for a generic perturbation 
the maximum growth is achieved by modes transverse to the magnetic field lines of the basic state. On the basis of 
this result, we have performed the multiple scale analysis for the marginally unstable case and for purely transverse 
perturbations. The analytic procedure yields a couple of partial differential equations which describe the full nonlinear 
evolution of magnetic perturbations. It is possible to show that these equations possess a Lyapunov functional and 
thus their solutions asymptotically approach a fixed point which represents a nonlinear equilibrium different from the 
basic one. The kinetic perturbations are always linearly stable, and their growth is uniquely due to the nonlinear 
coupling to the magnetic field. Numerical simulations of two-dimensional MHD performed with a pseudospectral code 
reveal an excellent quantitative agreement with the first-order analytical prediction in a wide range of values of the 
perturbative parameter. 

We therefore conclude that the loss of stability of parallel magnetic field configurations at moderate Reynolds 
number is due to the growth of large scale perturbations, and that the features of this instability can be captured by 
a multiple scale analysis. The transverse large-scale instability is likely to be the generic mechanism of instability of 
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sheared magnetic fields even for large Lundquist/Reynolds numbers whenever the basic state admits a large number 
of neutral lines. When there is a single neutral line at large enough Lundquist /Reynolds numbers this mechanism is 
overcome by the formation of resistive boundary layers. 
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FIG. 1. Marginal stability lines for magnetic flux and stream function. Region /: stable with T$ < and r$ < 0. Region 
77: r* > 0, r$ < 0. Region III: < 0, r 4> > 0. Region IV: instable with > 0, Rj, > 0. 
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FIG. 2. Growth rates 7 of the transverse Fourier modes k for simulation with r — 1/64, v = 0.49 and r\ — 1.0. The 
continuous line represent the linear prediction (|l9|). 
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FIG. 3. Time evolution of magnetic potential of the first Fourier transverse components of wavenumber k n = n/16 for the 
DNS with r = 1/16, ij = 0.4 and v = 1.0. The number of unstable modes is 4. 




O 2 71- O 2 7T 



FIG. 4. Snapshot of the magnetic flux tp (left) and stream function ip (right) for t = 200 (upper) and t = 20000 (lower). 
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FIG. 5. Time evolution of the first square Fourier components of solution of the Cahn-Hillard equation. Compare with 
Figure relative to the direct numerical simulation of MHD equations. 
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FIG. 6. Snapshot of the stream function tp reconstructed according to (Fit) for t = 200 (left) and t = 20000 (right). 
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